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Abstract 

We consider a type of Quantum Electro-Mechanical System, known as the shuttle system, first 
proposed by Gorelik et al. , [Phys. Rev. Lett., 80, 4526, (1998)]. We use a quantum master equation 
treatment and compare the semi-classical solution to a full quantum simulation to reveal the dynamics, 
followed by a discussion of the current noise of the system. The transition between tunnelling and 
shuttling regime can be measured directly in the spectrum of the noise. 
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I. INTRODUCTION 



Nanofabrication techniques, combined with single electronics, have recently enabled po- 
sition measurements on an electromechanical oscillator to approach the Heisenberg limit 
In this paper we present a master equation treatment of a version of a quantum electromechan- 
ical system (QEMS), the charge shuttle, first proposed by Gorelik 4 . In the original proposal a 
metallic grain is surrounded by elastic soft organic molecules and placed between two elec- 
trodes. This forms a Single Electron Transistor (SET) with a movable island. The coupling 
between the vibration of the island and the tunnelling onto the SET island dramatically alters 
the transport properties of the SET. The tunnelling amplitudes between the reservoirs and 
the island are an exponential function of the separation between island and the reservoirs. If 
the island is oscillating with a non negligible amplitude, this separation is a function of the 
displacement of the island from equilibrium and thus the tunneling current is modulated by 
the motion of the island. When there is a non-zero charge on the island the applied electric 
field accelerates the island. As the electron number on the island is a stochastic quantity, 
the resulting applied force is itself stochastic, but constant for a given electron occupancy 
of the island. Assuming the restoring force on the island can be approximated as harmonic, 
we have a picture of a system moving on multiple quadratic potential surfaces, with differ- 
ing equilibrium displacements, connected by conditional Poisson processes corresponding to 
tunneling of electrons on and off the island. The shuttle thus provides a fascinating exam- 
ple of a quantum stochastic system in which electron transport and vibrational motion are 
strongly coupled. 

In this paper we idealise the island to a single quantum dot with only one quasi-bound 
electronic state. This corresponds to an extreme Coulomb blockade regime in which the en- 
ergy required for double occupancy is not bound. This minimal model captures the essential 
quantum stochastic dynamics of the shuttle system. The quantum dot jumps between two 
quadratic potential surfaces, displaced from each other, corresponding to no electron on the 
island and one electron on the island. As noted by previous authors, the system exhibits 
rich dynamics including a fixed point to limit cycle bifurcation in which the average electron 
occupation number on the island exhibits a periodic square wave dependence. In this paper 
we give a quantum master equation treatment of this quantum stochastic dynamical system, 
with particular attention to the shuttling and the current noise spectrum. We use the Quantum 
Optics Toolbox 5 to compare and contrast the well known semiclassical predictions to the full 
quantum dynamics. In particular, we compare the picture of ensemble averaged dynamics of 
various moments with a 'quantum trajectory' 6 simulation of moments. A quantum trajectory 
is a concept taken from quantum optics to describe the conditional dynamics of the system 
conditioned on a particular history of stochastic events. Such conditional dynamics provide 
insight into the effect of quantum noise on the the semiclassical prediction of regular electron 
shuttling on the limit cycle. 

Various versions of a charge shuttle system have been experimentally investigated. A re- 
view of the theoretical and experimental achievements in shuttle transport can also be found 
in the work of Shekhter et al? ' . When a voltage bias is applied between the electrodes, a 
current quantisation resulting from electron interactions with the vibrational levels for dif- 
ferent voltage bias was found. By using C(,o embedded between two gold electrodes, Park 
et a/.— have demonstrated that indeed there is current quantisation for various bias voltage 
which results in a stair-like feature within the current-voltage curve. Although because of 
its high frequency (around Terra Hertz) and low amplitude oscillation, the molecule hardly 
shuttles between the electrodes in this setup, this experiment has provided key evidence of 
the involvement of vibrational levels in changing the properties of the current. This quan- 
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tized conductance also was observed in several other experiments 9- . Zhitenev et al. utilize 
metal single electron transistor attached on the tip of quartz rods as scanning probe while the 
experiment by Erbe et al.—, combines a nanomechanical resonator with an electron island to 
produce a QEMS system. The experimental setup used by Erbe is similar to the one proposed 
by Gorelik 4 . Huang et al. also reported the operation of a GHz mechanical oscillator 11 . 

Several attempts to explain the behaviour of the system have been offered both from clas- 
sical and quantum point of view. The current quantisation and its low frequency noise was 
investigated via a classical approach by Isacsson 12 . The current-voltage relation in the shuttle 
system exists within two regimes. The first regime is when the electron tunnels straight into 
the dot from the source and off to the drain, without much involvement of the island move- 
ment. This is called the tunnel regime. The C6o system lies within this tunnel regime. The 
other regime is when the island oscillates to accommodate the current flow, which we call 
the shuttle regime. However, measurement of average current alone cannot provide enough 
information to distinguish whether the system is in the shuttle regime or tunnel regime. It 
was shown that a calculation of the noise is needed in addition. Therefore the noise signature 
was first obtained by finding the Fano factor at zero frequency 13 . Recently Flindt et al. 14 have 
calculated the current noise spectrum using a method different form that used in this paper. 
We compare the two methods in section VI. 

Another interesting property of the system is the existence of a dynamical instability with 
limit-cycle behavior which was found in a similar setup using a single metallic grain placed 
on a cantilever between two electrodes^. This forms a three-terminal contact shuttle system. 
Classical analysis of the system points to the fact that this instability in the system leads to 
deterministic chaos. The semiclassical dynamics of the simpler case of the isolated island, 
the subject of this paper, was thoroughly investigated by Donarini et al.^* 

One of the early attempt to investigate the system within the quantum limit is given by 
Aji et al— where electronic-vibrational coupling is investigated both in elastic and inelas- 
tic electron transport by looking at the current-voltage relationship and conductance. Other 
properties of the transport within the shuttle system such as negative differential conductance 
have also been found 18 although the derivation only considers terms linear in the position of 
the island. Various conditions, such as when the electron tunnelling length is much greater 
than the amplitude of the zero point oscillations of the central island, have been investigated 
by Fedorets 19 . Using phase space methods in terms of Wigner function Novotny et al . 2 ^ 21 
identify crossover from tunnelling to shuttling regime. 

Another variation of the shuttle is offered by Armour and MacKinnon 22 . In this model 
the steady state current across a chain of three quantum dots system (one dot connected to 
each leads and one dot as vibrating island) was analysed by looking at the eigenspectrum. 
Numerical simulation here considers 25 phonon levels, within the large bias limit. 

In a recent thesis of Donarini 16 , the single dot quantum shuttle and the three dot shuttle 
system was investigated using Generalized Master Equation approach using Wigner distribu- 
tion functions. The current and Fano factor at zero frequency is also investigated. 

H. THE MODEL 

The system consists of a quantum dot 'island' moving between two electrodes, the source 
and the drain. This is analogous to a quantum dot SET in which the island of the SET is 
allowed to oscillate and thus modulate the tunnel conductance between itself and the reser- 
voirs. However unlike a SET we do not include a separate charging gate for the island. When 
a voltage bias is applied between the two electrodes, the electron from the source can tunnel 
onto the island and as the island moves closer to the drain the electron can tunnel off, thus 
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producing a current. Here we assume that only one electronic level is available within the 
island, a condition of strong coulomb blockade. 




FIG. 1: Schematic representation of shuttling between a source and a drain through a quantum dot. 

The electronic single quasi-bound state on the dot is described by Fermi annihilation and 
creation operators c,c', which satisfy the anti commutation relation cc^ + c^c = 1. While 
the vibrational degree of freedom is described by a displacement operator x which can be 
written in terms of annihilation and creation operators a and a^, with the commutation relation 
aa^ — a^a=\. 



a + a 1 ). (1) 



2mV 

The Hamiltonian of the system is given by: 

H =h® I c Jf c + U c n 2 (2) 

+ hva^a (3) 

+ h& s ka\a k + hto d kb\b k (4) 

— eE x c^c (5) 

+ £(r i ,, J E_(i) % c t + h.c)+£(r^ J B + (x)v t +h.c) (6) 

k k 

+ Y,s( at dp + adJ ? )+Y,h(o p dldp, (7) 
p p 

where E is the electric field seen by an electron on the dot. 

The first term of the Hamiltonian describes the energy of a single-electron quasi-bound 
state of the island. For the purpose of our simulation, we will scale other energies in terms of 
this island energy and thus conveniently set h&i = 1 . The Coulomb charge energy, U c is the 
energy that is required to add an electron when there is already one electron occupying the 
island (h = c' c). This energy is assumed to be large enough so that no more than one electron 
occupies the island at any time. This is the Coulomb blockade regime. In this regime it is 
better to regard the island as a single quantum dot rather than a metal island and we will refer 
to it as such in the remainder of this paper. The free Hamiltonian for the oscillator is described 
in term © where v is the frequency of the mechanical oscillation of the quantum dot. The 
electrostatic energy of electrons in the source (s) and drain (d) reservoirs is written as term 
©. With a l an d bfr, b, the annihilation and creation operator for the electron in the source 
and drain respectively. Term © describes the electrostatic coupling between the oscillator 
and charge while term © represents the source-island tunnel coupling and the drain-island 
tunnel coupling. In the shuttle system, the island of the SET is designed to move between the 
source and the drain terminal with an amplitude or fluctuation comparable to the distance of 
the island to the lead. Thus we introduce the term 

E ± (x) = e ±x/X (8) 

_ g ±T|(aW) (9) 
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with 

1/2 1 



^I.Sfj X (10) 

to account for the change in the tunnelling rate to the left and the right lead as the position of 
the shuttle varies. 

The last term, ©, describes the coupling between the oscillator and the thermo- 
mechanical bath responsible for damping and thermal noise in the mechanical system in the 
rotating wave approximation 23 . We include it in order to bound the motion under certain bias 
conditions. 

We now obtain a closed evolution for the system of quantum dot plus oscillator by tracing 
out over the degrees of freedom in the leads. A Markov master equation for the island- 
oscillator system can then be derived in the Born and Markov approximation using standard 
techniques 24 . If we assume the vibrational frequency of the oscillator is slow compared to 
bath relaxation time scales, we arrive at: 

p = — iv[a^a,p] 
+ jx[(a + a t )c t c,p] 

+ Yl(/(*(B/ -n L )<D[JE„{x)}p + (1 -f(h(D, - hl))(D [c£_(i)]p) 
+ J R (f(hm -ii R )(D[JE+(j£)]p + (1 -/(/i©/ -ii R ))<D[cE+(j£)]p) 
+ K(n p + l)'D[a]p + KHp2)[a ! }p, (11) 

with % = eEx\k and n p is the mean phonon number for the vibrational damping reservoir. We 
also have defined 

©[A]p=ApA t -^(A t Ap + pA t A), (12) 

where / is the Fermi function /(e) = l/(e e ' Td + 1). This Fermi function has an implicit 
dependence on the temperature, T e \, of the electronic system and the bias conditions between 
the source and the drain. The terms JliJr describe the rates of electron tunnelling form the 
source to the dot and dot to drain respectively. We have implicitly ignored co-tunnelling and 
higher order scattering events, so this equation applies under weak bias and weak tunnelling 
conditions. The final two terms proportional to K describe the damping of the oscillator, 
where h p = l/(e^ v /^ T — 1) and T are respectively the mean excitation and the effective 
temperature of a thermal bath responsible for this damping process. Thermal mechanical 
fluctuations in the metal contacts of the source and drain cause fluctuations in position of the 
center of the trapping potential confining the island, that is to say small, fluctuating linear 
forces act on the island. For a harmonic trap, this appears to the oscillator as a thermal bath. 
However such a mechanism is expected to be very weak. This fact, together with the very 
large frequency of the oscillator, justifies our use of the quantum optical master equation (as 
opposed to the Brownian motion master equation) to describe this source of dissipation 23 . 

In order to discuss the phenomenology of this system we first consider a special case. 
Under appropriate bias conditions and very low temperature, the quasi bound state on the 
island is well below the Fermi level in the source and well above the Fermi level in the drain. 
The master equation then takes the " zero temperature" form 

p = -rv[a t a,p]+i%[(a + a t )c t c,p] 
+ y L v[c*E-(x)]p+y R v[cE+(x)]p) 

+ K(n p + l)T)[a]p + Kn p <D[a*]p. (13) 
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The terms proportional to y^ and y^ describe two conditional Poisson processes, 
dNi{t),dN R {t), in which an electron tunnels on to or off the island. The average rate of 
these processes is given b y 2 V 6 i 27 

<E{dN L {t)) = y L Tr[E-(Jt)c f pcE-(jt)]dt, 
T,(dN R {t)) = y R Tr[E+(jt)cpc f E+(jt)]dt. 

where •£ refers to a classical stochastic average. Using the cyclic property of trace and the 
definition in Eq.© we see that 

<E{dN L {t)) = jLie-^cc^dt, (14) 
<E(dN R (t)) = y R (e 2 * lx Jc)dt. (15) 

It is now possible to see that the current through the dot will depend on the position of the 
oscillator. Under appropriate operating conditions (discussed below) we can use this depen- 
dance to configure the device as a position sensor or weak force detector. For a symmetric 
case where the tunnel-junction capacitances are almost the same, Q, ~ C R (neglecting the 
position dependence of the capacitances), the Ramo-Shockley theorem indicates that the av- 
erage current in the circuit can be given by 



/(f) = E (*(*)) = ! 



(16) 



If T| -C 1 and Jl = y R = y we may write this as 



where 



I(t) « e y/2 + ^(i(c t c-cc t )) + 0(f 2 ) (17) 
= ey/l + ^.-^ + ^x 2 ), (18) 



(x) k = Tr osc [x(k\p\k)} (19) 



with k = 0, 1 the occupation number states for the dot, and osc indicates a trace with respect 
to the oscillator Hilbert space alone. It is apparent that (x)k refers to the average position of 
the oscillator conditioned on a particular occupation of the dot. Clearly the average current 
through the system depends on the position of the oscillating dot. However the dependence 
on the first moment of position may be very weak. If the tunnel rates through the dot are 
much larger than all other time scales we expect that the occupation of the dot will reach an 
equilibrium value of \ quickly. In this case the term linear in position will be very small, 
leaving only a quadratic dependence. However if it can be arranged that yi ^ Jr, there will 
be a direct dependance of the current on the oscillator position. To clarify this situation we 
first look at a semiclassical description of the dynamics. 



IH. SEMICLASSICAL DYNAMICS 

The master equation Eq. (fTT1) enables us to calculate the coupled dynamics of the vibra- 
tional and electronic degrees of freedom. The equations of motion for the occupation number 
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on the dot and the average phonon number are 



d(c^c) 
dt 



+ Y*[/*«ccV i A> - (1 -f R )(c^ce^ x )}, (20) 



^ = ym 2 [f L (cSe-^) + (l-f L )(Sce-^)] 



+ Y*r| 2 [/*(ccV^> + {l-f R ){Jce- 2 *l 
-i%((a — a ! )c*c)+Kn — K(a ! a), (21) 

where the Fermi factors are defined by f a = /((O/ — /Jo) with a = L,R and /j a is the chemical 
potential in the source (a = L) and drain (a = R) and h(Di is the energy of the quasi bound 
state on the dot. The equation of motion for the average amplitude is relatively simple: 



d ( a ) ■ .i \ l / \ , • / t 



tv(a)-lK(a)+ix(c*c) (22) 



which is the equation of motion for a damped oscillator with time dependent driving. Unfor- 
tunately these first order number moments are coupled into higher order moments generating 
a hierarchy of coupled equations. A semiclassical approximation to the dynamics may be de- 
fined by factorising moments for electronic and vibrational degrees of freedom. This discards 
quantum correlations and thus is certainly not the appropriate way to describe a quantum lim- 
ited measurement. However it does enable us to see the essential features of the dynamical 
character of this problem. We will return to the full quantum problem in the next section. 

We begin the semi-classical approach by factoring moments of oscillator and electronic 
coordinates, for example of (cc*E_) into (cc^) (E 2 ), to obtain 

d -^ 1 = ld(E 2 )f L -^c)(E 2 )) 

+y R ((E 2 )f R -(Sc)(E 2 )). (23) 



Using the definitions, 



h 

x = T|A,(a-|-fl ) p = — i—, r(« — a ) 

2t|A 



we can write the semiclassical equations in terms of position x = (x), momentum p = (p) 
and electron number n = (c'c), 



jde- 2x/X f L - ne- 2x ' X ) + J R (e 2x / X f R - ne 2 *^) (24) 



dn 
dt 

dx p K 

17 = ~~^ x (25) 
dt m 2 

dp ? K /- 

— = —mvx p + Xy2mvhn (26) 

dt 2 

where we have made the further factorisation (E±(x)) = g ±2r A. These results agrees with 
the previous classical equations obtained by Isacsson 15 , in the case of zero gate voltage on 
the island. We will carefully consider the regime of validity of these semiclassical equations 
in section V. For now we note that factorising vibrational and electronic degrees of freedom 
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ignores any entanglement between these systems, while factorising the exponential assumes 
the oscillator is very well localised in position. 

In the zero temperature limit and appropriate bias we have that /l = 1, /r = 0. The 
semiclassical equations of motion then take the form 

§ = n^-n)e-^ x - lR ne^ x (27) 

d(X K 

— = -iva- -a + ixn (28) 
at 2 

with 

a=(a) = (x)/(2lr\) + i(p)lr[/h = X + iY . 

The system of equations, Eq . (12711281 has a fixed point, which undergoes a hopf bifurcation. 

To see this we begin by scaling the parameters by v and T|; ~ — > y, ~ — > K and — > % and 
v — > 1 by scaling time x = vt and redefining X and Y by letting a = r\(X + iY). Then 

7=TL(l-«K 4X -y^ (29) 
ax 

da k 

— — = —i(x — — oc + i%n. (30) 
ax 2 

The fixed point is given implicitly by 



y L e-^*+y R e^* 1 + ^e 



(3D 



i + (!) 2 



(33) 



from which we can see that it must satisfy, 



X = X,(1 + (^ 2 )(1 + ^). (34) 

At the hopf bifurcation the fixed point looses stability and a limitcycle is created. To see 
this, first obtain the linearized matrix about the stationary point. 

-A, 
DF = | -f 1 

X -i -I 



where 



The stability of the fixed point is determined by the eigenvalues of this matrix. If one 
or more of the eigenvalues have positive real part the fixed point is unstable. For complex 
eigenvalues the transition between stable and unstable occurs when the eigenvalues are pure 
imaginary. Here this is when 

_ A,k(A*(A, + k) + 1 + (| ) 2 ) 
8YlYa 
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At % = %h the eigenvalues are — (A* + k) , ±i/j where fj = J A*K + 1 + (^) 2 and the fixed point 
has a one dimensional stable manifold and a two dimensional center manifold. For % < %/, the 
fixed point is stable and for % > %/, it is unstable. This suggests a hopf bifurcation, however 
it is necessary to work out the stability coefficient to determine if it is subcritical, creating an 
unstable limitcycle or supercritical, creating a stable limitcycle. This involves some algebra. 
First transform the system in the vicinity of the fixed point to normal form via the matrix of 
eigenvectors P. 

-A*K(A* + f) -A*(A* + §) -A*(^ + l + (f) 2 ) 
Then in normal form coordinates u = P 1 [n — m*, X — X*, Y — Y*) T the system becomes 



du 

dx 



-(A*+k) 

—in | u + gNl (u) + SjlJ R ((% - % h ) (u\ + (35) 

o in o 



where g and f are column vectors, whose entries are gi = P n , fj = P a ■ Nl(u) is a scalar 
nonlinear function of w, obtained by perturbation. To cubic order in Uj 



M(u) = -4/i'X' y'A 2 - 4y L y* - Sn'X ,z A, ^ 
where 

(n / ,x / ,y / ) r = ^u T - 

Now the limitcycle bifurcates into the center manifold which is tangent to the u\ = plane. 
So if mi = h(u2i M 3) is the equation of the center manifold through (0, 0, 0) at % = %/„ then 
/j(0, 0) = and |^(0, 0) = 0. This means that a Taylor series approximation to the center 
manifold will have no constant or linear term and so the first nonzero terms are of quadratic 
order in ui and 

h(u2, ui) = a2Qu\ + a\\U2Uj, + ao2"3 + higher order terms, 

for some a2o, «n and ao2- Now differentiating u\ = h(u2, K3) gives; 

du\ dh du2 dh dui 
dx du2 dx dus dx 

On the center manifold ^(h(u2, W3), U2, W3) are functions of 112 and H3 only, so this 
equation can be used to calculate the coefficients ay in the Taylor series approx- 
imation to h(u2, U3) recursively, by equating coefficients of like powers of 112 and 
M3. Once h(u2,U3) is found this can be fed back into the equations of motion for 
U2 and uj, to obtain the approximate equations of motion on the center manifold. 
Finally the stability coefficient for a two dimensional system in normal form 2 ^ is 

a = if XXX ~\~ gxxy + /xyy + g>'yy) + JgQj ^ xx -^>'- v ) _ _ fxxgxx + fyySyy) 



evaluated at (0, 0), where here co = /j = yA*K + 1 + (| ) 2 . The subscripts indicate partial 
derivatives of function f or g with respect to the variables x and y. For instance fxxifxxx is 
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FIG. 2: Illustration of two possible type of Hopf bifurcation in the shuttle system with varying coupling 
X, a) supercritical and b) subcritical and saddle node bifurcation of the limit cycles. 
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FIG. 3: Plot of %/, for various fixed values of y# as a function of Yz,. The line is solid where the hopf 
bifurcation is supercritical and dashed where it is a subcritical hopf bifurcation. 



is a short hand for 2nd and 3rd derivative of function / with respect to x. Here the stability 
coefficient must be calculated numerically because the position of fixed point is only known 
implicitly via Eq. (l34b . 

Figure |3]plots %/, for various fixed values of as a function of Yl. The line is solid where 
the stability coefficient is negative, implying a supercritical hopf bifurcation and dashed 
where it is positive, implying a subcritical hopf bifurcation. 

At a supercritical bifurcation a stable limit cycle bifurcates from the fixed point, existing 
for %> %h- At a subcritical an unstable limitcycle bifurcates from the fixed point, existing 
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for %< Xh- Continuity of solutions as the parameter y L is changed suggests that the stable 
limitcycle existing for X> Xh above the solid critical line also exists above the dashed line. 
Numerical evidence shows this to be the case and that it continues to exist well below the 
dashed line, eventually being annihilated in a saddle-node bifurcation with the unstable limit 
cycle created in the subcritical hopf bifurcation at the dashed line. A schematic diagram of 
the two bifurcations are shown in figure El For Jr = 5 and ji = 1.01 the hopf bifurcation 
occurs at Xh = 0.42650636 and the saddle node bifurcation at % sn = 0.315. A glance at Fig. 
reffig:chih, where a vertical grey line indicates the range of % for Jr = 5 and y^ = 1.01 for 
which there are two limit cycles shows that there is a significant parameter region, where two 
limit cycles coexist. 

In general for fixed y^ the stability coefficient is positive for small and very large Jl and 
negative in between. This means that if ji and Jr are about 1, say, a stable limit cycle 
bifurcates and is present for x > Xh- But if (Jl/Jr) is much less than 1, a more complicated 
situation may arise for % < Xh, where an unstable limit cycle exists close to the critical point 
surrounded by a stable limit cycle. 

We then solved numerically the full system of equations, Eqs. (l27l) and d28t . for various 
values of the parameters. In the shuttling regime the electron number on the dot n[t) exhibits a 
square wave dependence as a single electron is carried from source to drain, where it tunnels 
onto the drain and the dot returns empty to the source to repeat the cycle. This is shown 
as the thin line in FigUta). The effect of shuttling generally occurs when the maximum 
displacement of the island is quite large, and where the strength of the tunneling depends 
strongly on the position of the island (k small). During shuttling, the electron number on the 
dot is constant. This gives, from Eq. d27t . an implicit relation between the shuttle position 
and the dot occupation, 

y L<? " 4T l* 

»W = ^ +Y ^- (36) 

Near the equilibrium point, X = 0, this implies that for Jl = Jr, n = 0.5. Away from the 
equilibrium point we have that 

»W={? 07) 

This behaviour is evident in the semiclassical dependance of n(t) (thin solid line) in FigEta). 

A condition for shuttling is given also by Gorelik 4 by specifying the requirement for the 
amplitude of the shuttle oscillation to be much bigger than the tunnelling length X. Donarini^ 
set the shuttling condition as to when the mechanical relaxation rate is much smaller than the 
mechanical frequency and also that the average injection and ejection rate is approximately 
equal to the mechanical frequency of the oscillator. 

The quantum dynamics may be determined by solving the master equation in the phonon 
number basis of the oscillator and the charge basis for the dot. It is necessary to truncate the 
phonon number basis high enough to include the amplitude of the limit cycle. 

To overcome the numerical difficulties with simulating large number of phonon levels for 
the quantum case described in Sec.V later, we choose a set of values of % and r\ which will 
give a rather small limit cycle in the semiclassical approximation in Fig|8] The accuracy of 
the semiclassical simulation is dependent on X as can be seen in Sec.V by comparing the 
factorized and unfactorized result from the numerical method. 

We now return to consider the dependance of the total current on the oscillator position. 
The total current through the device is given by Eq. (fT6l) . In the semiclassical approximation 
this is given by 

I T (t) = ^(l-n)e-^ x + J j-e^ x n. (38) 
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At the fixed point region, this n can be substituted by given in Eq. d3TT) to give: 



lj y L e- 4 nX*+y R e^ x *' (39) 
When T) is small, we can simplify the current further to: 

I T = 1^- (40) 

jL + yR-4r[X*(y L -y R ) 

Here we need to remember that the tunelling rates Jl and Jr determine the steady state posi- 
tion of X*. We can express, from Eqs. (I3TT) and d32l the tunneling rates as: 

lR X*(l+4rpY*) ' { > 



where for simplicity we have set: 



We can thus rewrite the current: 



S =?TW (42) 



h =^ B : x '\ («) 

B + 4V[X, y 

We can see that when T] is small the current lj is linearly dependent on the fixed point position 
X*, with a slope of 

We check this result using the full quantum simulation (Eqs. (l20"l) . (I2TT) ') and compare it with 
the result of the semiclassical current Eq. (l39b . We plot the result for various combination of 
T] and % in FigH For each condition, we vary the ratio of y# to Jl to give the plotted curve. 
As we can see from Fig. Stb),(c), the current is indeed linearly proportional to the steady 





FIG. 4: Current versus steady state position of the oscillator for various combination of r\ and % with 
varying ratio of Jl and Jr along each curve. Here we have chosen K = 0.2 a:r| = 0.3, X = 0.3, b: 
r| = 0.03, X = 0.3, c: r\ = 0.03, % = 0.5. Bold lines are results from the full quantum simulation and 
thin lines for the semiclassical approximations. 

state position of the oscillator when T| small. In this case, the semiclassical expression of the 
current given above in Eq. (l43t is a very good approximation of the actual current. 
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IV. A POSITION TRANSDUCER SCENARIO 



In this section, for simplicity, we assume that the zero temperature limit applies for which 
bound state of the dot is well below the Fermi level in the source and well above the Fermi 
level in the drain. The irreversible dynamics are then conveniently described in terms of 
two conditional Poisson jump processes with rates defined in Eqs. (ll4ll5l> . The jump process 
Eq. JLH) can only occur if there are no electrons on the dot, and the jump process Eq. (fT5]) can 
only occur if there is an electron on the dot. In the case that there is no electron on the dot, 
the quantum dot moves in a quadratic potential centered on the origin. In the case that there 
is an electron on the dot, the non-zero electrostatic force means the quantum dot oscillates 
in a quadratic potential displaced from the origin by Xq = %/v. We thus have a picture of 
a system moving on one or the other potential surfaces interrupted by jumps between them. 
This is schematically illustrated in Fig|5J Due to the exponential dependance of the jump 
rates on position (see EqsIT4landfT5l). the process dNi(t) is vastly more likely to occur when 
X < and conversely, the jump process dNi{(t) is much more likely to occur when X > 0. 
This means that the jump processes are an indication of which side of X = the dot is located. 



0\ n=1\^ 
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FIG. 5: A schematic illustration of the two potential surfaces connected by Poisson jumps. 



With this interpretation we can easily describe the conditional dynamics of the shuttle 
conditioned on a history of jump processes. In quantum optics such conditional dynamics are 
called quantum trajectories 6,29 . Let us suppose that at time t = t\, the occupation of the dot 
is zero and the jump process dNiif) occurs at t = t\ + dt\. The dot then becomes occupied 
while the state of the oscillator changes according to 2 ^ 

lYftfc)) ^ WL(tk + dt k )) = * c-^lvfa)) (44) 

\/pL[M,tk) 

where pLi/Jjk) = (V(//t)k ~ fa \ y V(tk)) an ^ we have defined /j = 2/X. With these definitions 
we see that £ (dN^t^)) = YLPL{n,tk)dt. We can develop some useful insight into what this 
state transformation means in the case that |\|/(^)) is a Gaussian with mean position of x& and 
variance (5k- In this Gaussian case we have 

where (2m) ! = 2. 4. 6.... 2m. After the jump process the mean position changes to 

{■yL\x\y L )=x k -2c k /X. (46) 



13 



This equation applies equally well to jumps to the right, dNR, with a change in the sign of X. 
Thus we see that if there is jump due to dNi, on average the conditional state moves to state 
with a mean closer to the source, while if a jump occurs to the right, cINr, the conditional 
state changes to a state with a mean position closer to the drain. This conditional behaviour is 
consistent with the interpretation of the jumps as effective measurements of the position of the 
quantum dot. More discussions on the quantum trajectory picture and numerical simulations 
on the conditional dynamics will be presented in the next section. 

V. SOLVING MASTER EQUATION NUMERICALLY 

With the help of the Quantum Optics toolbox^, we can solve the master equation directly 
by finding the time evolution of the density matrix. This was done by preparing the Liouvil- 
lian matrix in Matlab and solving the differential equation given the initial conditions. 

The expectation values for any desired quantities such as the electron number (c'c), the 
phonon number (a* a), position (x) and the momentum (p) of the oscillator can be calculated 
by tracing the product of this quantities with density matrix p. The result can then be plotted 
against time. The same method can be applied to calculate the steady state solution of the 
expectation values using p ss . 

The initial state of the system has been set up to incorporate the two electron levels, namely 
the occupied and empty state, combined with an N levels of phonon. The number of phonon 
levels included determines the accuracy of the calculation. Of course the more phonon levels 
included the more accurate the simulation will be. However only a limited number of phonon 
levels can be considered. This is due to the limited computer memory that is available and 
also considering the calculation time which will be significantly higher for larger N. Thus we 
try to find the minimum number of phonon levels which gives convergent results. This will 
ensure that the simulation still has a reasonably accurate solution. Donarini 16 use the Arnoldi 
iteration 30 to find the stationary solution of the matrix to overcome this memory problem. 
However here we have proceeded without, in the hope of looking at not only the stationary 
solution but also the dynamical evolution of the shuttle. 

The behaviour of the shuttle depends strongly on the rate of electron jump between the 
island and the leads. We investigate this by looking at the variation in the electron number 
expectation (c^c) at various rates JliJr- This is shown in Fig. 

reffig:3DVarGamma in which we have set Jl to be equal to Jr. When Jl,Jr are small, the 
electron number slowly increases until it reaches the steady state condition. In the region 
where the values of Jl^Jr is close to the frequency of the island, oscillation starts to occur, and 
depending on the damping that was set, the electron number can reach a steady oscillation 
putting the system well in the oscillatory regime. When Jl, Jr are very large compared 
to other frequency scales in the system, we will arrive at the strongly damped regime of 
the shuttle (see Sec. V.B), where the jump rate of the electron is fast enough to damp the 
oscillations in the electron occupation number of the island. Since we set to be equal to Jr 
the steady state happens at (c^c) = 0.5. 

Similarly the behaviour of the shuttle also changes according to T|, as described in Fig|7] 
At T] = the electron occupation number grows to a steady state. As we increased T] further, 
the oscillations start to occur with increasing amplitude. Here we use 100 phonon levels for 
the numerical calculation. 
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time(2jc) 

FIG. 6: Plot of the electron occupation number in the island for various tunneling rate Yl in a logarith- 
mic scale when r\ = 0.3 % = 0.5, v = 1 k = 0.05. 




time(2jt) 



FIG. 7: Plot of the electron occupation number in the island for various values of r\ when Yl = Jr = 1, 
X = 0.5, v = 1 k = 0.05. 

A. Oscillatory Regime 

The oscillatory regime will occur when the oscillation caused by the electron jump rate 
introduces continued kicks on the island. This happens when the jump rate is close to the 
oscillation frequency of the island (Yl ~ v). By setting an appropriate damping to this oscil- 
lation (k), there exists a condition where the island will keep oscillating between the leads. 
With this setup, the system will be in the shuttle regime. 

We then choose a set of parameters where the system shows the behaviour of a shuttle, 
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(a) (b) 



FIG. 8: QO toolbox simulation for shuttling condition (bold lines), compared with the semiclassical 
simulation (thin lines). Parameters are Yz, = Yr = 1,V = 1,T| = 0.3,% = 1,K = 0.05.. We set the initial 
condition as Xq = 4. 1 and Fo = to start the evolution close to the limit cycle, (a) Evolution of average 
values for position (solid lines) and momentum (dashed lines), electron number, current in the source 
and the energy of the oscillator; (b) Limit cycle behaviour in 3D. 

that is a continued oscillation of the electron number along with the oscillation of the island 
position. To ensure the convergence of the numerical solution, we use a smaller value of T| 
that will still give a shuttling behaviour. We choose a combination of T] and % that will give 
the smallest limit cycle to minimise the truncation error. 

Within the region where the limit cycle exists, we can plot the electron expectation num- 
ber against the average position and momentum and observe the shape of the limit cycle. We 
explored both the full density matrix simulation and the semiclassical solution to be com- 
pared. From the result (FigEl), the quantum simulation appears to be more damped than its 
semiclassical counterpart. This is due to the effect of the noise. This slight difference can 
also be caused by the dependence of the electron number on its correlation with the position 
that was ignored in the semiclassical case. 

To check this we have plotted the difference between the factorized and unfactorized mo- 
ment at this particular variable combination (Fig©. The time range in which the difference 
in the factorized and unfactorized occurs agrees well with the time range when the semiclas- 
sical and the quantum simulation disagree in Fig. 

reffig:qo shuttle. This disagreement happens at the time when the shuttle is in transition be- 
tween the zero and one electron occupation number. 

Of course the truncation will pose some inaccuracy in the quantum simulation at a longer 
time. However we have checked that this is not the case at least for a short period of time by 
comparing it with a simulation that includes a larger phonon number. 

To investigate the effect of T| on the correlation between the factorized and unfactorized 
moments, we can plot (cc^a + a^)) and (cc^fa + cfl). We can see from FigO that the 
semiclassical approximation agrees with the quantum simulation under the condition that T| 
is small enough. As r\ increases, the evidence of this difference becomes noticeable. This 
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FIG. 10: Plot of the difference between the factorized and unfactorized moments for various r\ with 
Yl = Y« = 1,X = 0.5,k = 0.05. 

difference oscillates with peaks located at times when electron jumps happen, that is when 
the oscillator is near the equilibrium position. 

As opposed to what the name "electron shuttle" suggests, the dynamics in the shuttle 
regime, for the parameters specified in Fig 03 is not like a conventional shuttle which picks 
up an electron when it is closest to the source and drops the electron when it is closest to 
the drain, as also suggested by Nord et alQ. Looking at the rate of the average electron 
number and the average current in the source (FigEta)), this is certainly not the case. The 
shuttle picks up an electron near an average displacement of zero, slightly towards the source, 
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and continues to travel closer to the source electrodes. It then oscillates back and drops the 
electron at a slightly displaced average position from equilibrium towards the drain electrode. 
The shuttle then continues to get closer to the drain before oscillating back to repeat the cycle. 

An important distinction must be made between the dynamics of the averages derived from 
solving the master equation and a dynamics conditioned on a particular history of tunneling 
events. This distinction is already suggested by inspecting the average electron number as a 
function of time. In any actual realisation of the stochastic process, the number of electrons 
on the dot is either zero or one, yet the ensemble average occupation number varies smoothly 
between zero and one. The reason for this is that the actual times at which transitions between 
the two states takes place fluctuates. 

We can more easily appreciate this distinction using an alternative approach to understand- 
ing the dynamics based on 'quantum trajectories'. The quantum trajectory method (some- 
times called the Monte Carlo method) first introduced in quantum optics, is a method of 
looking at the evolution of a system conditioned on the results of measurements made on that 
system. 2 ^ 2 ^ 27 ' 2 ?'? 2 . This method will allow one to monitor 'events' such as the jump of an 
electron to the island which causes the displacement kick on an oscillator. 

The Quantum Optics Toolbox enables a direct computation of the conditional dynamics of 
the operator moments by implementing a so called 'jump unravelling' of the master equation. 
First we plot a sample trajectory for a slow electron jump rate Jl — 0.1 to see the effect of 
electron jump on the evolution picture of the system. A random jump of electrons from the 
source to the island (FiglTTI) according to rate Jl.Jr was introduced. The dynamics of the 
shuttle as a position transducer, as predicted in section IV, can be seen in the conditional 
averages of the displacement. The variable r\ controls the amount of displacement of the 
island when an electron jumps on and off onto the island. Larger value of T| caused a larger 
displacement kick when a jump occurs. During the time when the electron is on the island, 
the phonon number of the oscillator oscillates with a similar behaviour to the oscillation of 
the position. 

The single trajectory for the shuttle case with the same parameter in FiglUcan be seen 
in FigO The electrons mostly jump onto the island from the source when it is closer to 
the source and jump off when closer to the drain. At the jump, the island gets a slight dis- 
placement kick towards the source when jumping on and towards the drain when jumping 
off. However this does not stop the shuttling motion of the island and does not repel it to 
the opposite direction as suggested earlier by Nord et al.^>. It can also be seen that when the 
electron manages to jump onto the island when island is still close to the drain, it is more 
probable for the electron to jump off straight away. 

The conditional dynamics of the system just described corresponds to an experiment in 
which number of electrons on island is monitored continuously in time. As we can see, the 
behaviour of the conditional dynamics differs from the behaviour of the ensemble average. 
However, averaging over many different realization of the trajectories as shown in FiglTZl 
would lead to a closer and closer approximation of the ensemble average behaviour in FigOD 

B. Strongly Damped Regime 

There are two ways of damping the shuttle into the fixed point regime. One is to damp 
the motion of the shuttle itself by introducing a large mechanical damping K. Alternatively 
we can damp the oscillation of the electron occupation number in the island. This happens 
when the rates of the electron jump Jl,Yr are large compared to the natural frequency of the 
island vibration. The fast electron jumps act as an internal damping to the shuttle. Within this 
regime the electron number expectation (c^c) monotonically approaches 0.5 when ji = Jr. 
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FIG. 11: Plot of a single trajectory showing the dynamics of the jump and the result in the phonon 
number. Parameters are y L = y R = 0.1, V = 1 ,r| = 0.3, X = 1,K = 0.05. 
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FIG. 12: Plot of a single trajectory showing the dynamics of the jump when Yl = Yr = 1,v = 1,T| = 
0.3, % = 1,k = 0.05. 
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When the bare electron tunnelling rates are very large compared to other frequency scales 
in the problem, we may assume the dot approaches its steady state for bare tunnelling quickly 
as compared to the typical time scale of the oscillator. In this case, p(t) = p s ^ x p (t). The 
bare can be substituted into the density matrix of the total master equation and then be 
traced out with respect to the dot degrees of freedom to get an effective master equation which 
involves only the reduced density matrix of the oscillator. This effective master equation can 
be calculated from the reduced density matrix, from Eq. (fTTTl . 

Since r\ is assumed to be small, we can expand the expression to second order in T|: 

e^ x = 1 + r\X + (r\X) 2 / (2!) H . We can then re-write the zero temperature full master 

equation as: 

p = — iv[a^a,p] 
+ 2ix[Xn,p] 
+ 2y L (l-n)r| 2 [X,p£,p]] 
+ 2y R nT\ 2 [X,[X,p]] 

+ K(n p + l)®[a]p + Kn p 'D[a ( }p, (47) 
with n = Yl/Yl + Y^- The following moments can thus be derived from Eq. d47t : 

d(a^a) 



2xn(Y) -4Yl(1 - n)T} 2 - 4y R nT} 1 + Kn - a) (48) 



dt 



^=v(F)-|(l) (49) 
d{f) v<*> (50) 



dt 2 

The moments (a) and (Y) form a closed system of differential equation which can readily 
be solved 



(jt) = e- {K/2)t \^Xo-X*)cos(vt) + (Yo-Y*)sm{vt)J+X*, (51) 
(?) = e - {K/2)t ( -(X -X*)sm(vt) + {Y -Y*)cos(vt) ) +Y*. (52) 



where again X* and Y* is simply the displacement in the equilibrium such as given in Eqs. d32l) 
and d33l) with = n. Xq and Yq is the initial condition of X and Y respectively. The analytic 
expressions of Eqs. (l5Tb and d52l are useful for checking the solution of the master equation 
given by Matlab, to ensure that the truncation in the phonon number is adequate. 

The shuttle oscillation is damped to the new displaced position of X* which agrees to the 
obtained result previously. When Yl = Y# = gamma, we have = n = j in the regime when 
the tunneling rates are very large compared to other frequency scales (especially when T| is 
relatively small). In this case, the oscillatory behaviour of Eqs. (BTT) and (l52l) do not depend 
on the actual values of the tunneling rates gamma. It can also be deduced that the decay rate 
of the oscillation envelope is e _K? / 2 . In this regime, the result of the analytical expressions 
matches the quantum simulation quite well. 



C. Co-existence regime 

As discussed in Sec. Ill, we can also have a regime in which the behaviour of the shuttle 
depends on its initial condition. The system will either be attracted to the limit cycle and thus 
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FIG. 13: Average evolution of the shuttle with a large Yi and y#. Here we have chosen Yz, = Yr = 
10, r) = 0.1, % = 0.1, K = 0.05,Xo = 0.5, Fo = for both quantum simulation with N = 25 (thick lines) 
and semiclassical solution (thin lines). 

be in the shuttle regime or be attracted to the fixed point and be in the tunneling regime de- 
pending on its initial condition within the correct parameters where the subcritical bifurcation 
occurs. Following previous authors 16 , we call this the 'co-existence regime'. 

Semiclassically this can be seen when we plot the average evolution of the shuttle. De- 
pending on the initial conditions, the shuttle will either be attracted to the fixed point position 
or undergoes the stable limit cycle oscillation (Fig. [141) . 




FIG. 14: Semiclassical limit cycle of the shuttle in the co-existence regime. Here we have chosen 
Yl = 0.03, Yr = 1,11 = °-3, X = 1, K = 0.05, X Q = 2.2 and Y = 0. The evolution starts close to the 
unstable limit cycle and then moves toward the stable limit cycle 

The quantum average calculation in this regime however does not show the subcritical 
bifurcation since averaging over the noise in the system dampens this effect. This can be 
seen in the evolution of the single trajectory which is captured to the fixed point position at 
random times. A sample of trajectories each from different initial conditions were plotted in 
figure [15] and [T6J 
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FIG. 15: Average evolution of the shuttle. Here we have chosen Yz, = 0.01, y« = 1,T| = 0.3, X = 1,K = 
0.05, X = 1.2, F =0. 
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FIG. 16: Average evolution of the shuttle. Here we have chosen = 0.01, y/j = 1,T| = 0.3, X = 1,K = 
0.05, X = 5,F = 0. 

D. Finite Temperature 



We can easily extend these calculations to the finite temperature case by including the 
fermi factor fi and fy, which was previously set to 1 and for the zero temperature case, in 
the calculation for both the full quantum simulation and the semiclassical approximation. 

The effect on temperature on the system is shown in Fig. [T7] Comparing this with pre- 
vious result for the zero temperature (Fig. |8j there is a suppression of the electron number 
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FIG. 17: Average evolution of the shuttle with finite Temperature Yz, and y«. Here we have chosen 
j L = y R = i,n = 0.3, x = 1,k = 0.05, X = 4.1,1b = 0,/l = 0.8 and f R = at CO/ = 3v for both 
quantum simulation with N = 41 (thick lines) and semiclassical approximation (thin lines). 

oscillation. There is also a significant difference between the quantum and the semiclassical 
simulation for the electron number occupation which resulted in a difference in the current in 
each of the leads. 

VI. NOISE CALCULATION 

In surface gated 2DEG structures some recent experiments monitor the charging state 
of the dot via conductance in a quantum point contact 33 . However such techniques cannot 
easily be adapted for a nanoelectromechanical system. In experiments involving tunneling 
through a double barrier quantum dot structure the simplest thing to measure is the source- 
drain current. In the QEMS experiments of Park et al 8 and also of Erbe et al, the source 
drain current carried signatures of the vibration of the nanoelectromechanical component. In 
this section we calculate, using the Quantum Optics Toolbox, the current noise spectrum and 
show that it indicates the transition between the fixed point and the shuttling regime. 

The current seen in the external circuit, when electrons tunnel on and off the dot, only 
indirectly reflects the quantum nature of the tunneling process. Tunneling causes a local de- 
parture from equilibrium in the source and drain reservoirs that is restored through a fast 
irreversible process in which small increments of charge are exchanged with the external cir- 
cuit. While tunneling obviously involves a change of charge in units of ±e, the increments of 
charge drawn by the external circuit are continuous quantities determined by the overall ca- 
pacitance and resistance of the circuit. The current responds as a classical stochastic process 
conditioned on the quantum stochastic processes involved in the tunneling. In many ways 
this is analogous to the response of a photo electron detector to photons. 

The connection between the quantum stochastic process of tunneling and the current ob- 
served in the external circuit is given by the Ramo-Shockley theorem and is a linear combi- 
nation of the two Poisson processes defined in Eqs. <l 1 4-1 15b . The noise spectrum of such a 
current involves moments of both the tunneling processes, and correlations between them. In 
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Sun et al. , one can find a detailed example of how such correlations are determined by the 
corresponding master equation for the quantum dot system. 

Recently Flindt et al. 14 , have calculated a noise spectrum for the shuttle system defined in 
terms of the fluctuating electron number accumulating in the drain reservoir. Here we adopt 
a different (but equivalent) approach based on the framework of quantum trajectories. In this 
section we calculate, using quantum trajectory methods, the stationary current noise spectrum 
in the source current alone as this suffices to illustrate how the current noise spectrum reflects 
the transition from fixed point to shuttling. The total current shows the same features but has 
a different noise background. 

The two time correlation function quantifies the fluctuations in the observed current and 
is defined by: 

G(x) = e -i^(x)+E(I(t)I(t + x))]ti, 

The first term is responsible for shot noise in the current, while the second term quantifies 
noise correlations. We now show how the second term can be defined in terms of the station- 
ary state of the quantum dot itself. 

Let p(t) be the density operator representing the dot at time t. What is the conditional 
probability that, given an electron tunnels onto the dot from the drain between t and t + dt, 
another similar tunneling event takes place a time x later (with no regard for what tunneling 
events have occurred in the mean time)? If an electron tunnels onto the dot from the drain at 
time t, the conditional state of the dot (unnormalised), conditioned on this event is given by 

p(!)(t) =Y Le - f /Vp(?)c<r i/ \ (53) 
Given this state, the probability that another tunneling event takes place a time x later is 

G(t,x) = ^(e-^cJe^Xt)]) 

= yltr(e- M /hc f e m [e-^p{t)ce- £ / K ]\ . 

where formally we have represented the irreversible dynamics from time Mo t + x as the 
propagator e Lx . Let us now assume that the first conditioning event takes place at a time t 
long after any information about the initial state of the quantum dot has decayed away. That 
is to say the first conditioning event occurs when the dot has settled into the stationary state, 
Poo = lim f _>oop(0. The stationary two-time correlation function for the source current is then 
defined by 

G(x) = T^tr (V^cV^/VpooCe^]) (54) 

In terms of the dimensionless position operator, X, the noise in the two time correlation 
functions becomes 

G(x) =E(I L (t)I L (t + x))1>l =y z L TT[e- 4 ^cc^e L \e- 2 ^c^ 00 ce- 2 ^)] 

where e L% is the master equation evolution. 

The noise power spectrum of the current is given by: 

/»oo 

S((0)=2 G(x)(e ioyi + e- iaz )dx (55) 
Jo 

This noise spectrum can be directly calculated using the Quantum Optics Toolbox by first 
calculating the steady state solution poo and setting e~ 2T|X c^pooCe~ 2r|X as an initial condition 
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for the master equation evolution. Then we can calculate the expectation value of the operator 
cc^e~ Ar[X in the state evolved, according to the master equation, from this initial condition. 
It is important to note that the master equation does indeed have a steady state even in that 
parameter regime in which the semiclassical dynamics would imply a limit cycle. This is 
because quantum fluctuations cause a kind of phase diffusion around the limit cycle. These 
quantum fluctuations are precisely the random switchings observed in the single quantum 
trajectory shown in figure [121 hi fact as shown in 20 the Wigner function of the steady state 
has support on the entire limit cycle. 

The example of the noise spectra for various r\ is shown in FigHU The transition between 
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FIG. 18: Plot of spectrum in the left junction with various T|, setting % = 0.5. 

the tunnelling regime and the shuttling regime is clearly evident in the noise power spectrum, 
figure fT~8l (We have subtracted off the shot noise background). 

Setting T) = 0, % = we arrived back at the known spectrum for the source current in a 
double barrier device 34 with a single dip at zero frequency. As the values of T] and % (or ji) 
are increased, the frequency spectra develop sidebands which correspond to the frequency of 
the oscillator (Fig. 

ref fig : corrspectrum JT9l> . As the system approaches the shuttling regime, the frequency spectra 
pick up noise peak at zero frequency and additional peaks at higher frequencies close to a 
multiple of the oscillator frequency. This is a signature of the limit cycle formation. On the 
limit cycle, the frequency is shifted from the base oscillation frequency (0 = 1.055v. This 
is also given by the imaginary part of the eigenvalues of the linearised matrix expressed 
/a = a/A*k+ 1 + (k/2) 2 . This observation agrees with the predicted slight re-normalization 
of the frequency by Flindt et al.—. 

A similar feature of the noise is found by Armou^ in a system consisting of a SET that is 
coupled to a nanomechanical resonator. Although this is a different system from the shuttle 
system, the classical spectrum noise in this system also shows the dependency of the current 
on the position of the nanomechanical resonator. 
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FIG. 19: Plot of spectrum in the left junction with various Jl, setting r\ = 0.3, X = 0.5, K = 0.05. 
VII. CONCLUSION 

The dynamics of the shuttle system has been investigated via both the semi-classical and 
the full quantum master equation treatment. The latter reveals subtle properties of the dy- 
namics which was not found using the semiclassical treatment. The master equation is solved 
numerically using the Quantum Optics Toolbox enabling a detailed comparison of the semi- 
classical dynamics with the quantum ensemble averages. For the first time in the study of the 
quantum shuttle we compute the moments for the quantum state conditioned on a particular 
history of tunnelling events. This is called a quantum trajectory and it reflects what can be 
observed experimentally by monitoring the electron on the island. 

The conditional dynamics differs from the behaviour of the ensemble average, and gives 
new insight into the shuttling dynamics. In the shuttling regime, the ensemble average dy- 
namics of the electron occupation number is a smoothed square wave that slowly decays to a 
steady state value of one half. Given that the occupation number of the dot is either zero or 
unity this ensemble averaged behaviour may seem unexpected. However looking at the occu- 
pation number in a single conditional state (see Fig. fTTb indicates what is going on. A single 
quantum trajectory shows that the average occupation number is indeed either zero or unity 
and in the shuttling regime behaves like a square wave for short times but, at random times, 
suffers a phase jump. The ensemble average of many such trajectories with phase jumps at 
random times leads to the observed ensemble average dynamics as computed from the master 
equation. These random phase jumps ultimately lead to a steady state density operator for 
the system that, in the Wigner representation, is diffused around the limit cycle, as noted by 
Novotny et al.— . 

The shuttle dynamics was investigated in two regimes: the fixed point and the shuttle 
regime. In the fixed point regime, the shuttle is damped to a new displaced position. We have 
shown that there is a strong relation between the current and the fixed point of the position. 
This relationship is linear when the tunnel length is large (r| small). Thus it is possible to 
use the shuttle in a position transducer scenario. In this regime, the semiclassical treatment is 
shown to be accurately sufficient to describe the dynamics. 
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We provide the condition in which the shuttle regime will appear from the system by iden- 
tifying the appearance of limit cycle in the phase space of the shuttle. A careful analysis of 
the nonlinear dynamics using centre manifold method indicates that when Jl = Jr, the limit 
cycle forms through a supercritical pitchfork bifurcation. However when ji ^ Jr there is a re- 
gion of parameter space in which the bifurcation can be subcritical, and for which hysteresis 
is possible. Adjusting the damping K with respect to these parameters will cause the shuttle to 
be sufficiently damped and thus allow the shuttling to take place. The shuttle regime also ap- 
pears when the rate of the electron tunnelling is close to the oscillator frequency. The shuttle 
regime corresponds to the continuous oscillation of the electron number and results in addi- 
tional peaks at multiples of the limit cycle frequency in the noise spectra. This is destroyed 
when K is too large or when a large electron jump y^, Jr are introduced to the system. Both of 
these conditions will damp the shuttle into the displaced equilibrium position. The quantum 
shuttle thus provides a fascinating example of a quantum stochastic system in which electron 
transport is coupled to mechanical motion. In future studies we will investigate how such a 
system can be configured for sensitive force detection. 
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